Unstable attractors induce perpetual synchronization and desynchronization 
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Common experience suggests that attracting invariant sets in nonlinear dynamical systems are 
generally stable. Contrary to this intuition, we present a dynamical system, a network of pulse- 
coupled oscillators, in which unstable attractors arise naturally. From random initial conditions, 
groups of synchronized oscillators (clusters) are formed that send pulses alternately, resulting in a 
periodic dynamics of the network. Under the influence of arbitrarily weak noise, this synchroniza- 
tion is followed by a desynchronization of clusters, a phenomenon induced by attractors that are 
unstable. Perpetual synchronization and desynchronization lead to a switching among attractors. 
This is explained by the geometrical fact, that these unstable attractors are surrounded by basins 
of attraction of other attractors, whereas the full measure of their own basin is located remote from 
the attractor. Unstable attractors do not only exist in these systems, but moreover dominate the 
dynamics for large networks and a wide range of parameters. 



As attractors determine the long-term be- 
havior of dissipative dynamical systems, the 
notion of attractors is central to studies 
in many fields of science. According to 
the mathematical definitions of an attractor, 
states that originate in a certain volume of 
state space, the basin of attraction, evolve to- 
wards the respective attractor. Since states 
slightly perturbed from an attractor often 
stay confined to its vicinity and finally re- 
turn to it, attractors are widely considered to 
be stable. Attracting yet unstable states are 
consistent with an attractor definition intro- 
duced by Milnor. There is evidence that such 
Milnor attractors might not be uncommon in 
certain systems that exhibit strange invariant 
sets with a fractal geometry. In general, how- 
ever, unstable attractors seem to be special 
cases that have to be constructed artificially 
by precisely tuning parameters. Contrary to 
this intuition, we present an example of a dy- 
namical system, a network of pulse-coupled 
oscillators, in which unstable attractors with 
periodic dynamics arise naturally. Unstable 
attractors are shown to prevail for large net- 
works and a wide range of parameters. In the 
presence of arbitrarily weak noise, a perpet- 
ual synchronization and desynchronization of 
groups of oscillators occurs which leads to a 
switching among different attractors. 



I. INTRODUCTION 

The concept of attractors is underlying the analy- 
sis of many natural systems as well as the design of 
artificial systems. For instance, the computational 
capabilities of neural networks are controlled by the 
attractors of their collective dynamics ||l|, |[ ^. Con- 



sequently, the nature and design of attractors in such 
systems constitute a focus of current research inter- 
est g 1 1, 0, |, I 0, |ll In gen- 
eral, the state space of a nonlinear dynamical system 
is partitioned into various basins of attraction from 
which states evolve towards the respective attractors. 
Since states that are slightly perturbed from an at- 
tractor often stay confined to its vicinity and eventu- 
ally return to the attractor, attractors are commonly 
considered to be stable [TgI p7| p|. 

In this paper we study the dynamics of networks 
of pulse-coupled oscillators g ^, ^, in which un- 
stable attractors exist and arise naturally as a col- 
lective phenomenon. Such models of pulse-coupled 
oscillators describe e.g. synchronization in spiking 
neural networks and the dynamics of other natural 
systems as diverse as pacemaker cells in the heart, 
populations of flashing flreflies, and earthquakes (cf. 

61 1, 1, §, 0, 0, @ 0, g 111, 0, a 13, 11). 

We identify an analytically tractable network ex- 
hibiting unstable attractors. For this network we 
demonstrate the existence of attractors that are un- 
stable and located remote from the volume of their 
own basins of attraction. Such attracting yet unsta- 
ble states are consistent with a definition of attrac- 
tors introduced by Milnor, which neither presumes 
nor implies stability |Q. In certain other systems 
such Milnor attractors might not be uncommon if 
these systems exhibit strange invariant sets with a 
fractal geometry |6|, |2^. More gener- 

ally, however, attractors that are not stable seem 
to be special cases that have to be constructed ar- 
tificially by precisely tuning parameters. Contrary 
to this intuition, in the system considered here, un- 
stable attractors with regular, periodic dynamics are 
typical in large networks and persist even if the phys- 
ical model parameters are varied substantially. The 
first discovery of the occurrence and prevalence of 
unstable attractors in networks of pulse-coupled os- 
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dilators was reported in [y^. Here we give a more 
detailed analysis of such unstable attractors and ex- 
plain the observation that they only occur for exci- 
tatory (phase advancing) interactions but are absent 
if the interactions are inhibitory (phase retarding). 



We argue that dynamical consequences of unstable 
attractors may persist in a general class of systems 
of pulse-coupled units. Such consequences include 
an ongoing switching among unstable attractors in 
the presence of noise. In systems where the con- 
vergence towards an attractor has a functional role, 
such as the solution of a computational task by a 
neural network ||l|, ||, |^ , switching induces a high de- 
gree of flexibility that may provide the system with a 
unique advantage compared to multistable systems: 
In general, it is hard to leave a stable attractor after 
convergence, e.g. the completion of a task. With an 
unstable attractor, however, a small perturbation is 
sufficient to leave the attractor and to switch towards 
another one. 



This paper is organized as follows: In section || we 
introduce a class of network models of pulse-coupled 
oscillators. We briefly review earlier results on syn- 
chronization phenomena in such networks within the 
framework introduced by Mirollo and Strogatz that 
is described in detail in section III. In section |^ 
we describe the numerical observation that, in the 
presence of noise, trajectories approach and retreat 
from periodic orbits by perpetual synchronization 
and desynchronization of groups of oscillators. To- 
gether with further numerical investigations, this 
leads to the hypothesis, that the observed dynamics 
is induced by attractors that are unstable. In section 
we perform an exact stability analysis of a partic- 
ularly selected set of attractors. It is followed by an 
analysis of the dynamics during a switching transi- 
tion between two states (Sec. VI) that is further cor- 
roborated by a numerical investigation of the struc- 
ture of basins of attractions in state space. These 
results demonstrate that unstable attractors indeed 
exist in the class of networks of pulse-coupled oscil- 
lators considered. Section VII completes our analy- 



sis showing that unstable attractors prevail in large 
networks for a wide range of parameters. The paper 
concludes in section VIII with (i) a discussion of the 
dynamical consequence of switching among unstable 
attractors in comparison to smoothly coupled sys- 
tems, (ii) a brief presentation of preliminary results 
for networks exhibiting different interactions or more 
complex structures and (iii) an outlook for future in- 
vestigations. 



II. MODELS OF PULSE-COUPLED 
OSCILLATORS 

In studies of synchronization phenomena in net- 
works of pulse-coupled oscillators, single elements 
are often modelled as phase oscillators, assuming 
that the dynamics of the amplitude of the oscillation 
is less important. Thus there is only one relevant dy- 
namical variable, Wi, that describes the dynamics of 
an individual oscillator i. In models of many natu- 
ral systems, the oscillator variable Wi represents an 
analogue of a potential, like the membrane poten- 
tial in the case of a current-driven nerve cell. The 
dynamics of a network of pulse-coupled oscillators is 
commonly described by a system of coupled ordinary 
differential equations 



dt 



A{W^) + B{W^)S,{t) 



(1) 



for i G {l,...,Af} where A and B are continuous 
functions and Si represents the interactions within 
the network. The pulse-coupHng between oscillators 
is given by 

N oo 



where Sij is the strength of the coupling from oscil- 
lator j to oscillator i and the response kernels Kij (t) 
have the property that Kij{t) > 0, Kij{t) = for 
t < 0, and J^^Kij{t)dt = 1. The sum includes all 
times tj^m at which oscillator j reaches a threshold 
We := 1 (from below) for the m^^ time. 



W,{t,^m) > Wg = 1 



dW,{t) 



dt 



> 0. (3) 



If this threshold is reached, Wj is reset to a potential 



W,{t 



+ ) 



w,. 







(4) 



and a signal is generated that is sent to oscillators i. 
Commonly, a unit that is reset and sends out a pulse 
when it reaches a threshold is said to "flre" at that 
instant of time. This form of pulse-coupHng ideal- 
izes the fact that in diverse biological systems such as 
populations of flashing flreflies or networks of spiking 
neurons in the brain, units interact by stereotyped 
short-lasting signals that are generated as the state 
of a unit reaches a threshold. Note that the normal- 
izations Wg = 1 and Wreset = are made without 
loss of generality. Whereas Eqs. (|l])-(Q) describe the 
interaction dynamics without delays between send- 
ing and reception of a pulse, a delay r > can easily 
be included by a transformation Kij (t) Kij {t — r) 
of the response kernels. 
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It is often convenient ||8| (and under weak condi- 
tions possible) to transform the variables Wi accord- 
ing to 

V^{t) ~ - y B{w)-^dw (5) 

where 



C= / B{w)-^dw 

10 



such that (ih becomes 



(6) 



(7) 



for i £ {l,...,iV} where A{V) 
C-^A{W{V))/B{W{V)) is determined solving 
Eq. (|) for W. If the duration of the response of a 
unit i to an incoming pulse is sufficiently brief, the 
kernels Kij in (||) may be idealized by the Dirac 



distribution 



K,,{t) = S{t) 



(8) 



such that the coupling becomes discontinuous. 

For the commonly used leaky integrate-and-fire 
models of oscillators the free dynamics is given by 
the Hnear inhomogeneous differential equation 



dV 
dt 



-jV + I 



(9) 



such that the network dynamics of such oscillators is 
given by substituting AiVi) = I — jVi in (|^ where 
/ is an external current and 7 > measures the dis- 
sipation in the system. This linear differential equa- 
tion has the advantages that the solution is known 
explicitely and the response to an additional cur- 
rent can be found by superposition arguments. For 
sufficiently large external current, / > 7, the free 
{Si{t) = 0) dynamics 



V{t) 



7/7(1 
V{t 



-7M 



T) 



for < t < T 
Vit) 



(10) 



is periodic (Fig. |^) with period T=iln(l — ^) ^ . 

To study the synchronization of pacemaker cells 
in the heart, Peskin studied a simple, globally cou- 
pled network of such leaky integrate-and-fire oscilla- 
tors ||l9| where the coupHng is not delayed, r = 0, 
and the responses are infinitely fast, Kij{t) — S{t), 
excitatory, and homogeneous, = e > 0. In his 
1975 book he conjectured that arbitrary initial con- 
ditions converge towards the fully synchronous state 
in which all oscillators fire simultaneously. He gave 
a proof for N — 2 oscillators assuming that both 




Figure 1: Dynamics of a non- interacting integrate-and- 
fire oscillator. Without threshold, the potential V would 
converge to its long-time limit I /y (dashed line). Every 
time the threshold Ve = 1 is reached, the potential V is 
reset to zero and a pulse is sent. 



coupling strength and dissipation are small, e ^ 1, 
7 < 1. 

In a seminal work of 1990, Mirollo and Strogatz 
generalized the approach of Peskin . Contrary to 
Peskin, whose analysis was based on the Hnearity 
of the differential equations, Mirollo and Strogatz's 
only assumptions were that the free {Si{t) = 0) dy- 
namics can be described by a phase-like variable and 
that the interactions are mediated by some potential 
function U that is a monotonically increasing and 
concave down function of this phase |^ (for model 
details see below). This framework includes many 
cases for which the associated differential equation 
is nonlinear such that its solution may not be known 
explicitely and superposition arguments fail. Within 
this general framework, they proved that, for all TV, 
almost all initial conditions will ultimately end up 
in the fully synchronous state. For systems without 
dissipation (7 — 0), equivalent to Hnear functions 
U, Senn and Urbanczik ||lj] proved that the dynam- 
ics becomes fully synchronous even if the intrinsic 
frequencies and the thresholds of the oscillators are 
not quite identical. Like the previous investigators, 
they treated the case without delay, r = 0, for which 
the fully synchronous state stays the only attractor. 
This periodic orbit is an example of period-one dy- 
namics, for which every oscillator fires exactly once 
during one period, and is the simplest dynamics such 
a system may exhibit. 

Within the framework introduced by Mirollo and 
Strogatz, it was, moreover, shown that the introduc- 
tion of a delay time t > 0, that occurs ubiquitously 
in natural systems, changes this situation drastically 
|§[ ^: With increasing network size, an exponen- 
tially increasing number of attractors coexist. In a 
large region of parameter space, these are periodic 
orbits with period-one dynamics, that exhibit several 
groups of synchronized oscillators (clusters), which 
reach threshold and send pulses alternately. It was 
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also shown ^ that these kinds of attractors arise 
not only in excitatorily coupled networks (e > 0) but 
also if the oscillators are coupled inhibitorily (e < 0) . 

In the following part of this paper, we show that 
for excitatory coupHng, many of these periodic or- 
bits with period-one dynamics are unstable attrac- 
tors. As a consequence, trajectories converge to- 
wards these unstable attractors by synchronization 
of groups of oscillators into several clusters but, in 
the presence of arbitrarily weak noise, diverge sub- 
sequently via desynchronization of clusters. 



III. MIROLLO-STROGATZ MODEL 

We consider a homogeneous network of N all-to-all 
pulse-coupled oscillators with delayed interactions. 
A phase-like variable ipiit) G (— oo, 1] specifies the 
state of each oscillator i at time t such that the dif- 
ference between the phases of two oscillators quanti- 
fies their degree of synchrony, with identical phases 
for completely synchronous oscillators. The free dy- 
namics of oscillator i is given by 



d(j)i/dt = 1. 
Whenever oscillator i reaches a threshold 

= 1 

the phase is reset to zero 

4>^{t+) = 



(12) 



(13) 



and a pulse is sent to all other oscillators j ^ i, which 
receive this signal after a delay time t. The interac- 
tions are mediated by a function U{(j)) specifying a 
'potential' of an oscillator at phase 4>. The function 
U is twice continuously differentiable, monotonically 
increasing, U' > 0, concave (down), U" < 0, and 
normalized such that U{0) = and U{1) = 1. 

For a general U((j)) we define the transfer function 



i?e(0) = C/"'(C/(0) 



(14) 



that represents the response of an oscillator at phase 
to an incoming subthreshold pulse of strength i 
that induces an immediate phase jump to 0"*" = 
Hi{if)). Depending on whether the input e is sub- 
threshold, U{4>) -I- e < 1, or suprathreshold, U{(f)) -\- 
e > 1, the pulse sent at time t (Eq. (|l2|)) induces a 
phase jump after a delay time r at time t' = t + t 
according to 



H,{^,{t')) if l7(,^,(t')) + £< 1 
if U{(l)j {t')) + i > 1 
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Figure 2: An incoming pulse of strength i induces a phase 
jump <P+ ~ 01 (t+) = U-\UiMt)) + ^) = Hi{,p^) that 
depends on the state 0i := </'i(t) of the oscillator at time 
t of pulse reception. Due to the monotonicity of U, an ex- 
citatory pulse (e > 0) induces an advancing phase jump. 
If the incoming pulse puts the potential above threshold 
{U{4>i) -I- £ > 1), the phase is reset to zero {(pf — 0). An 
inhibitory pulse (e < 0) would induce a regressing phase 
jump such that the phase may assume negative values 
(not shown). 



This phase jump (Fig. g) depends on the phase 
(l)j{t + t) of the receiving oscillator at a time r af- 
ter the signal by oscillator i has been sent at time t, 
(11) the effective coupling strength — e = s/{N — 1), 
and the nonlinear potential U. The interactions are 
either excitatory (e > 0) inducing advancing phase 
jumps (Fig. H) or inhibitory {e < 0) inducing re- 
tarding phase jumps. Note that in response to the 
reception of an inhibitory pulse, the phases of the 
oscillators may also assume negative values. 

The general class of functions U captures the dy- 
namics of a variety of systems. In particular, given 
any differential equation of the form (|^ the free 
{Si{t) = 0) dynamics of which has a monotonic pe- 
riodic solution V{t) with period T and negative cur- 
vature, the function U can be taken as the scaled 
solution. 



U{cl>) V{^T). 



(16) 



By this transformation, the general class of pulse- 
coupled oscillators (Eqs. (0)-(Q)) with infinitely fast 
response (^ is mapped onto a normalized phase de- 
scription (Eqs. (pl])-(p!5|)). For instance, the stan- 
dard leaky integrate-and-fire oscillator (||) with so- 
lution (fol) leads to 



UiAd^) - -(1 

7 



(17) 

Another example is given by the conductance-based 
threshold model of a neuron |^^, in which A(VF) = 
7(VFeq - W) and B{W) = g(Ws - W) with equilib- 
rium potential Wcq > 1, membrane time constant 
7 > 0, synaptic reversal potential Ws, and conduc- 
(15) tivity g > 0. After a transformation of variables 
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V = V{W) according to §) and a scaling we 
obtain 



In (1 - wr') 



(18) 

The resulting potential function U is always mono- 
tonically increasing and, for a wide range of biologi- 
cally reasonable values of Ws and Wcq, also concave 
down. 

For all numerical studies and simulations pre- 
sented in this paper, we use the functional form 

= t/fc(0)= In (l + (e^ -!)(/.) (19) 

that results in an afHne transfer function Hi{4>) = 
e^^cj) + const, that was utilized also in the original 
work of Mirollo and Strogatz [|| as well as in the 
later investigation of the influence of delay ^ . 

Compared to systems of nonlinear differential 
equations, the Mirollo-Strogatz approach has the ad- 
ditional advantage, that numerical calculations can 
be performed exactly on an event-by-event basis. 
Given the state of the system at time i, defined by a 
phase vector 4){t) and a list of signals sent together 
with their sending times, the dynamics can be com- 
puted numerically by iterating between two kinds of 
events (in an appropriate order): 

1. If the next event is the reception of a signal 
after a time Ati, shift all phases 4>i{t) by this 
amount and apply the map according to 
Eq. ( p^ . If the received signal is subthreshold 
for oscillator j, its resulting phase is given by 



(l)j(t + Mi)^ Hi{(f)j{t) + Ati). 



(20) 



if the signal is suprathreshold for oscillator 
j' , 4'j'{t) + A^i > 1, its phase is reset to zero. 



0y(i-l- Ati) = 0.0, 



(21) 



and a signal is generated, i.e. its sending time 
t + Ail and the set of sending oscillators j' is 
stored. If A: > 1 signals have been simulta- 
neously sent by synchronized oscillators their 
simultaneous arrival can simply be numeri- 
cally realized by an enlarged coupling strength 
{k — l)e or fee, depending on whether or not 
the considered receiving oscillator has sent one 
of the k signals. 

2. If the next event is that one or more oscillators 
reach threshold after time Ai2 — 1— max^ 4>i{t) 
(without receiving an additional input signal 
after a time Ati < Ai2), those oscillators j are 
reset to zero. 



and a signal is generated (see event 1.) whereas 
the phases of all other oscillators j' are just 
translated in time according to 



,/.jv(t + At2) = 0j'(i) + At2. 



(23) 



^j(t + Ata) =0.0, 



(22) 



These are the only two possible events that change 
the state of the system, such that the simulation time 
is increased by Ati or Ai2, depending on the kind of 
event: To determine the state at the next time step 
after time t numerically, one first finds the minimum 
of the time Aii after which the next signal would 
arrive and the time At2 after which the next phase 
would cross threshold without an additional incom- 
ing signal (cf. Eq. ([l^)). Thus, time steps smaller 
than AiMS = min{Aii , At2} do not need to be con- 
sidered. Note that AtMS depends on the state of the 
system such that it changes with time t. 

This event-based numerics can be exploited for ev- 
ery choice of [/((/>), and is not restricted to C/if(0) 
derived from the linear differential equation (^, for 
which one may as well exactly numerically integrate 
the original differential equation due to linear super- 
position (cf. BI^). 

Hence, in contrast to standard numerical integra- 
tion of nonlinear differential equations, for which the 
time steps A^de have to be taken sufficiently small 
and the numerics is only approximate, the above nu- 
merical algorithm for the Mirollo-Strogatz model is 
exact. It is also fast, if oscillators constitute syn- 
chronized groups, because then typical time steps 
are long, AiMS > A^de- 



IV. PERPETUAL SYNCHRONIZATION 
AND DESYNCHRONIZATION 

For such pulse-coupled systems, periodic orbits 
with groups of synchronized units constitute rele- 
vantattractors g |, | 0, |, |, 0, 0, |l|, |l|, |l|, 
poj , |2l| , P2| . As mentioned above, for the networks of 
Mirollo-Strogatz oscillators with delayed interactions 
(r > 0) considered here, many different cluster-state 
attractors with several synchronized groups of oscil- 
lators (clusters) coexist H, ||. After an initial tran- 
sient, these networks settle down onto such a periodic 
orbit that displays period-one dynamics with clus- 
ters firing successively within each period. Thus, in 
particular, phase-differences between oscillators are 
constant at all times when a reference oscillator is 
reset. 

Whereas the noise-free dynamics is seemingly simi- 
lar for both kinds of coupling, in the presence of noise 
the collective behavior of excitatorily coupled oscil- 
lators strongly differs from that of oscillators cou- 
pled inhibitorily. For inhibitory coupling all cluster- 
state attractors are stable against small perturba- 
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Figure 3: Phase dynamics of a large network (A*' = 100, r — 0.15). Phases of all oscillators are plotted whenever 
a reference oscillator has been reset such that the axis labeled 'Time' is discrete and nonlinear. (a,b,c) Inhibitory 
coupling (e = —0.2). (d,e,f) Excitatory coupling (e = 0.2). (a,d) Dynamics with noise (noise level r] = 10~^). (b,e) 
Deterministic dynamics in response to a single phase perturbation (arrow, perturbation strength a — 10"'^); (b) 
for inhibitory coupling the system returns to the attractor; (e) for excitatory coupling the system switches from a 
six-cluster to a five-cluster state. (c,f) Phase differences from the average phase of one cluster (shown in b and e, 
respectively) in response to the perturbation. 



tions. Under the influence of sufflciently weak noise 
the system stays near some periodic orbit that has 
been reached after a transient from a random ini- 
tial state (Fig. ||a). The dynamics after a small per- 
turbation in an otherwise noiseless system conflrms 
this stability property: Perturbations to all clusters 
decay exponentially and the original attractor is ap- 
proached again (Fig. P3,c). 

In contradistinction, for excitatory coupling, we 
flnd (cf. Fig. ||d) that, although the system con- 
verges towards a periodic orbit from random initial 
states, small noise is often sufficient to drive the 
system away from that attractor such that succes- 
sive switching towards different attractors occurs. 
In principle, this dynamics might be due to stable 
attractors located close to the boundaries of their 
basins of attraction, such that noise drives the tra- 
jectory into a neighboring basin. If this explanation 
were correct, ever smaller perturbations off the at- 
tractor would lead to an ever lower probability of 
leaving its basin. In an otherwise noiseless system 
we tested this possibility by applying instantaneous, 
uniformly distributed, independent random pertur- 
bations 5i e [0, cr] of gradually decreasing strengths 
(down to tr = 10^^"*) to the phases 4>i of all oscillators 
i after the system had settled down to an attractor 
(see Fig. ||e,f). Even for the weakest perturbations 
applied, none of the perturbed states returned to 
the attractor, but all trajectories separated from the 
original attractor. We thus hypothesized that the 



persistent switching dynamics (Fig. |3|i) is due to at- 
tractors that are unstable. 



V. STABILITY ANALYSIS 

In order to verify this hypothesis directly, we an- 
alyze a small network of iV = 6 excitatorily cou- 
pled oscillators for which instantaneous perturba- 
tions lead to a similar switching among attractors. 
We flx parameters to e = 0.2, t = 0.15 and the po- 
tential function ?7(0) = C/3(0) according to Eq. (p^. 
At these parameters the network exhibits a set of 
periodic orbits with period-one dynamics that are 
related by a permutation of phases in such a way, 
that the system may switch among them (cf. Fig. ||a, 
points on the periodic orbits marked in red, yel- 
low, blue). These orbits are structurally stable in 
a neighborhood of these parameters and the func- 
tion U . Each attractor can be characterized by a 
list of cluster occupation numbers giving the num- 
ber of oscillators Uk in the fc'^ cluster, counted in 
the order of increasing phases at times after a refer- 
ence oscillator {i — 1) has been reset. For instance, 
the attractor marked in yellow is characterized by 
the occupation list [2,2,1,1] such that the transi- 
tions among the attractors (Fig. ^) are described by 
the sequence ([1,2, 2, 1] (red) -> [2, 2, 1, 1] (yellow) -> 
[2, 1, 1, 2] (blue) [1, 2, 2, 1] (red)) for the particular 
set of perturbations applied. For each of these three 



attractors there exists another permutation-related 
attractor with the same occupation Hst but with the 
phase values of the two clusters containing only one 
oscillator interchanged. It turns out that, depending 
on the perturbation, transitions from one attractor 
occur to one of only two other attractors that are 
related through the latter permutation. Thus, once 
the system has settled down to one of the attractors 
shown, it may switch within a set of six periodic orbit 
attractors in response to sufficiently small perturba- 
tions. 

Due to their permutation-equivalence these or- 
bits have identical stability properties. The state 
of the network at time t is specified by cj){t) = 
(01 (t), . . . , 6e{t)y , such that the orbit marked in yel- 
low in Fig. 0a is defined by the initial condition 



0(0) = (0, 0, A, A, B, C)T 



where 



A = H,{t) 

B = H2e{l + 2T-ai) 

C = H2e{Hi{2T) + 1 ^ 



and recursively defined 

a, = U-\he + U{T + ai^i)) 



(24) 



(25) 
(26) 
(27) 



(28) 



for i G {1, . . . ,4}, ao = 0, and ki — — k/^ ^ 1, 
k2 = 2. Here the origin of time, t — 0, was chosen 
such that oscillators 1 and 2 have just sent a signal 
and have been reset. Moreover, at i = only these 
two signals (and no others) have been sent but not 
yet received. The numerical values for the particular 
parameters considered, A w 0.176, B « 0.499, C w 
0.747, can be identified in Fig. ^ (orbit marked in 
yellow). This orbit indeed is periodic. 



</.(T) = </,(0), 



(29) 



and, in particular, period-one, such that each oscil- 
lator fires exactly once during one period T. 

To perform a stability analysis, we define a 
Poincare map by choosing oscillator z = 1 as a refer- 
ence: Let 



(30) 



be the perturbed phases of the oscillators i at times 
tn > 0, n e N, just after the resets of oscillator 1, 

0i(4)^O. (31) 

Thus the five-dimensional vector 

Sr.^<l)r.- (0, A, A, B, C'Y (32) 
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Figure 4: Small network (AT = 6, e = 0.2, r = 0.15) 
exhibiting unstable attractors. (a) Noise-free phase dy- 
namics in response to instantaneous perturbations of 
magnitude cr = 10~^ (arrows). When an attractor is 
reached, the phase configuration specifying the current 
cluster state attractor is marked in color. The pertur- 
bations induce a split-up of clusters and a divergence 
from the attractor such that the network reaches differ- 
ent attractors successively. For the realization of per- 
turbations shown, the attractors marked in red (cluster 
occupation list [1,2,2,1]), yellow ([2,2,1,1]), and blue 
([2,1,1,2]) are visited cyclically. For each of the three 
attractors there exists another attractor with the same 
cluster occupation list but with phase values of the two 
clusters with only one oscillator interchanged, (b) Basin 
structure of the three attractors color-marked in (a) in 
two-dimensional planar section through six-dimensional 
state space. The planar section is defined by one point on 
each of the three periodic orbit attractors, represented by 
small red, yellow, and blue disks, respectively. Basins of 
the three attractors are marked in the same colors. The 
part of state space within this section is of hexagonal 
shape. Darker gray areas are basins of the three other 
permutation-related attractors (cf. (a)) that are located 
outside the section shown. Lightest gray marks the union 
of the basins of all other attractors. 
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defines the perturbations Sn,i for i £ {2, . . . , 6} where 
we choose < Sn,2 and <5„_3 < (5„.4. Here 

C' = H2i{Hi{2T + Hi{Q)) + 1 + r - 04) (33) 

that numerically yields C" ~ 0.756. Because after a 
general perturbation clusters are split up such that 
in particular < 02, oscillator 6 receives the (later) 
signal from oscillator 1 only after it is reset by the 
(earlier) signal from oscillator 2 resulting in C ^ C. 
Thus, the original orbit is super-unstable, i.e. has 
infinite expansion rate and the stability analysis is 
performed for the orbit obtained in the limit ^„ ^ 0. 
It is important to note that unstable attractors also 
exist that do not posses such a "partner orbit" but 
have a finite expansion rate themselves. For conve- 
nience, we here continue considering the above set 
of orbits that allow a straightforward explanation of 
the phenomenon of unstable attraction. Following 
the dynamics, the five-dimensional Poincare map is 
given by 



where F is defined by 



(34) 



A 



F2{S) = 

F3(<5) = Li-L^ + Hiir-Li + Ls) 
F4{6) = Hg{T + L4 - L3) - A 
Fsid) = He{Hg{l + 2t - + L4 - L^) ~ B 
Fe{6) = Hg{L4 - L-^ + Hg{T + 1 - L3+ 

(35) 

with the abbreviations 

L, := L,{8) = Hs{t + Hg{T -82 + HS2 + 
Hs{T + S,~S2 + He{T))))) 

(36) 

for i e {3, 4}. Here the difference L^ — L^ in Eqs. ( |35| ) 
is of order L4 - L3 ^ 0{\\S\\). 
The linearized dynamics 



(37) 



of a slightly perturbed state with split-up clusters is 
described by the Jacobian matrix 



M = 



dF{S) 



dS 



s=o 



/O \ 

a -a ' 

0-13 f3 

* * * 

V * * * 00/ 



(38) 



where a,/3 > and * denote nonzero real numbers. 
It has four zero eigenvalues 



A, = for i G {1,2,3,4} 



(39) 



which imply that a six-dimensional state-space ac- 
cessed by the random perturbation is contracted 



Figure 5: Simultaneous suprathreshold excitatory input 
synchronizes immediately due to the reset at thresh- 
old. Sufficiently small but positive phase differences 
\<p2{t) — (pi(t)\ > are reduced to zero such that 
(ji^if^) = 4>i{t^) ~ 0. This is the mechanism of di- 
mensional reduction of effective state space, by which 
attractors can be built. 



onto a two-dimensional manifold. This reflects the 
fact that suprathreshold input received simultane- 
ously by two or more oscillators (Fig. ||) leads to 
a simultaneous reset and thus a synchronization of 
these oscillators independent of their precise phases. 
If a single oscillator is reset by a suprathreshold sig- 
nal, it instantaneously exhibits a precise lag in firing 
time At = T compared to the oscillator that has sent 
this signal. 



In particular, the zero eigenvalues (Eq. (|39|)) re- 
fiect the following contracting dynamics: (i) Per- 
turbations of phases 6n,5 ^ or 6n,6 ^ are re- 
stored immediately by suprathreshold input pulses 
received from oscillators j = 6 or j = 1, respec- 
tively. This gives rise to the eigenvalues Ai = A2 = 
corresponding to the eigenvectors vi oc (0, 0, 0, 1, 0)^ 
and V2 oc (0,0,0,0,1)^. (ii) A splitting of the first 
cluster, Sn,2 > = (5„,i , corresponding to the vector 
V3 oc (1, 0, 0, 0, 0)^, is restored after one period due to 
one suprathreshold input pulse from oscillator j = 4. 
At the same time, however, this splitting induces a 
perturbation of oscillators i = 5 and i — 6 such that 
<^n-i-i,5 7^ and Sn+i,e / which is restored accord- 
ing to (i) in the subsequent period. Taken together, 
this accounts for the eigenvalue A3 = 0. (iii) As long 
as (S„^3 = SnA , a perturbation vector is mapped onto 
the subspace spanned by vi and V2 (see (i)) within 
one period and is then mapped onto zero during the 
next period. As a result, the eigenvalue A4 = cor- 
responds to the direction V4 oc (0, 1, 1, 0, 0)^. In ad- 
dition to this analysis, a stability analysis for the 
subset of states with the cluster (j)3{t) = (f>4{t) kept 
synchronized, results in superstable directions only, 
as expected. 

In contradistinction to this contracting dynamics, 
the concavity of U implies that simultaneous sub- 
threshold input to two or more oscillators leads to 
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Figure 6: Simultaneous subthreshold excitatory in- 
put desynchronizes due to the concavity of U . Small 
phase differences \(l>2(t) — > are increased, 

\4>2{t'^) — 4>i{'t^)\ > \4>2{t) — 01 (i) I, providing the mech- 
anism that creates an instability. 

an increase of their phase differences, i.e. a desyn- 
chronization of oscillators with similar phases (Fig. 

For the orbits considered here, this is reflected 
by the only non-zero eigenvalue 

, ^ (2j7^(co) - U'{a,))U'{c,)U'{c2)U'ics) ^ 
' U'iai)U'ia2)U'ia3)U'{ai) 

(40) 

where 

Ci=T + ai (41) 

for i 6 {0, 1, 2, 3} and the are defined in Eq. (|28|). 
Because Ci > ai > Ci-i for all i and U' > 0, U" < 0, 
this eigenvalue is larger than one, i.e. the periodic or- 
bit is linearly unstable. This eigenvalue corresponds 
to a split-up of the cluster composed of the oscilla- 
tors i = 3 and i — A. Because the Jacobian (|3^ ) 
is not symmetric, the eigenvectors are not orthog- 
onal such that the corresponding eigenvector is not 
V cx (0, 1, —1, 0, 0)^ but has a component in this di- 
rection. If there is no homoclinic connection, this 
instability implies that such an attractor is not sur- 
rounded by a positive volume of its own basin of at- 
traction, but is located at a distance from it: Thus, 
every random perturbation to such an attractor state 
- no matter how small - leads to a switching towards 
a different attractor. 



VI. UNSTABLE ATTRACTORS 

Furthermore, this periodic orbit indeed is an at- 
tractor: According to the stability analysis, after 
two firings of the reference oscillator, a trajectory 
perturbed off a periodic orbit (e.g. the one marked 
in red in Fig. |^a, which is permutation-equivalent to 
the yellow one) is mapped onto a two-dimensional 
manifold, re-synchronizing one cluster. The trajec- 
tory then evolves towards a neighborhood of another 



attractor (here: the yellow one) in a lower dimen- 
sional effective state space without further dimen- 
sional reduction. Here, forming the second cluster, 
suprathreshold input leads to the last dimensional 
reduction while the state is mapped directly onto 
the periodic orbit. 

In general, a periodic orbit is unstable if, af- 
ter a random perturbation into its vicinity, one or 
more clusters are not re-synchronized by simultane- 
ous suprathreshold input but desynchronize due to 
simultaneous subthreshold input. An unstable at- 
tractor results if these clusters are formed through 
synchronization in a region of state space that is lo- 
cated remote from the periodic orbit towards which 
the state then converges. 

Although this is a discontinuous system with de- 
layed interactions such that there is no simple basin 
structure in the state space of phases and signals, 
a three-dimensional cartoon of the basin structure 
in a state space of phases may help to gain further 
insight about how trajectories approach and retreat 
from an unstable attractor in the presence of noise. 
Figure ^ shows that the basin volume is contracted 
by creating (at least) one cluster in a region of state 
space that is remote from the attractor itself. In con- 
trast, near the attractor, the same cluster is unstable 
against a split-up of the phases of the oscillators it 
contains. Basically, such an unstable attractor might 
be viewed as an unstable periodic orbit with a remote 
basin attached to its stable manifold that ensures the 
attractivity property. 

It is important to note that for inhibitory coupling 
we observe that all attractors are stable: An intuitive 
explanation is that there is only a mechanism of syn- 
chronization (Fig. ^ due to the concavity of U that 
contracts state space volume such that all attractors 
with period-one dynamics are stable. It is instructive 
to compare Fig. ^ for inhibition to Fig. || for excita- 
tion that display how simultaneously incoming sub- 
threshold pulses affect phase differences. Since for 
inhibition all attractors are period-one states |^, ||, 
we expect that the possibility of unstable attractors 
is excluded for this kind of coupling. 

In order to further clarify the structure of state 
space of networks of excitatorily coupled oscilla- 
tors, we numerically determined the basins of at- 
traction of the three attractors displayed in Fig. ^a 
in two-dimensional sections of state space. The ex- 
ample shown in Fig. ^ reveals that attractors are 
surrounded by basins of attraction of other attrac- 
tors as predicted by the above analysis. Because 
of this basin structure, noise induces repeated at- 
tractor switching among unstable attractors. Start- 
ing from the orbit defined by (|2j) the system may 
switch within sets of only six periodic orbit attrac- 
tors as is apparent from the basins shown in Fig. 0b. 
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cluster formation (synchronization) by suprathreshold input: discontiuous jump 



Figure 7: Cartoon of the basin structure of an unstable attractor A. First, a positive basin volume B{A) of states 
is mapped onto a lower-dimensional effective state space. This is achieved by simultaneous suprathreshold input to 
a group of oscillators, that synchronizes them to form one cluster. Because the attractor is located remote from its 
own basin volume the same cluster may desynchronize in response to small perturbations near the attractor, where 
incoming pulses are no longer supra- but subthreshold and thus lead to a desynchronization. 




Figure 8: Simultaneous inhibitory input synchronizes 
due to the concavity of U. Contrary to networks of os- 
cillators coupled excitatorily, here the inhibitory interac- 
tions allow subthreshold input only and lead to a decrease 
of phase differences, \(j)2{t^) — (l)i{t^)\ > \4>2{t) — 0i(i)| • 
Hence, these networks possess a mechanism for synchro- 
nization, but there is no simple possibility of desynchro- 
nizing a cluster state. 



However, in larger networks (of. e.g. Fig. ^) a clus- 
ter may split up in a combinatorial number of ways 
and exponentially many periodic orbit attractors are 
present among which the system may switch. The 
larger such networks are, the higher the flexibility 
they exhibit in visiting different attractors and ex- 
ploring state space. 

Until now, the analysis has focussed on a small 
network of = 6 oscillators, for which certain pe- 
riodic orbits have been demonstrated to be unsta- 
ble attractors. To study the desynchronization of 



clusters also observed in larger network in greater 
detail, we numerically determined the divergence of 
small random perturbations to an attractor in a net- 
work of A'' = 100 oscillators. As an example we 
chose two perturbations — aid* and ^2 = ^2^* 
into the same random direction 6* £ [0, 1]^ where 
(Ti = 10~^^ and 0-2 = (Ti-|-10~^'*. Figure ^ shows that 
the separation A :— maxj — (52,i| between the two 
perturbed trajectories exponentially increases with 
time. This indicates that also for large networks, 
desynchronization is due to a linear instability. Let 
us remark that, due to the spiltting-up of clusters 
by a general perturbation (cf. Eqs. (32) and (^3|)), 
two perturbations into independent directions might 
first lead to an additional discontinuous separation, 
followed by an exponential expansion. 

As a second quantity, that characterizes the dy- 
namics of switching, we consider the time needed 
by the system to switch from an unstable attrac- 
tor towards a different attractor after a random per- 
turbation of magnitude a. The perturbations ap- 
plied to the attractor are random phase vectors, 
drawn from a uniform distribution on [0, cr]^, with 
a e [10-12,10-2]. As displayed in Fig. |l^, for suffi- 
ciently small a, this switching time clearly increases 
exponentially with decreasing a. In particular, this 
indicates that the attractors found are indeed unsta- 
ble and do not possess small contracting open neigh- 
borhoods. Furthermore, it confirms our observations 
that the time of switching is mainly determined by 
the time of divergence from the original unstable at- 
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Figure 9: Separation of two perturbations off an unsta- 
ble attractor into the same random direction for a large 
network (iV = 100, e = 0.2, r = 0.15). The separation 
A grows exponentially with the time after the perturba- 
tion, measured in number of firing events of a reference 
oscillator. 
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Figure 10: Time of switching between two attractors de- 
pending on the perturbation strength a for a large net- 
work (iV = 100, £ = 0.2, r = 0.15). The switching 
time increases exponentially with decreasing perturba- 
tion strength. The discrete time axis (number of firing 
events of a reference oscillator) leads to a regular step- 
ping that slightly disrupts the exponential trend. 



tractor. Thus, in the presence of external noise, we 
expect a similar monotonic increase of an approxi- 
mate switching time with the amplitude of the noise. 

Taken together, this indicates that there exist also 
unstable attractors in larger networks that are en- 
closed by basins of other attractors. Thus, these 
results strongly support the hypothesis, that the 
switching found in large networks in the presence of 
noise (cf. Fig. H), is also due to unstable attractors. 
It is important to note that, without a perturba- 
tion, numerical noise does not induce a divergence 
of trajectories from attractors which are unstable: 
Synchronization occurs by simultaneously resetting 
the phases of two or more oscillators to zero (cf. Fig. 
H). Due to the global homogeneous coupling, all sig- 
nals simultaneously received by these oscillators are 
of the same size and are exactly synchronous nu- 
merically (see Sec. III). Thus, although the phase- 
advance computed for such signals is influenced by 



numerical round-off errors, such errors will be iden- 
tical for phases of synchronized oscillators and hence 
not induce a numerical desynchronization. 



VII. PREVALENCE AND PERSISTENCE 

The preceeding analysis demonstrates the exis- 
tence of unstable attractors. For excitatory cou- 
pling, these unstable attractors coexist with stable 
attractors. For a cluster-state attractor to be stable, 
all clusters necessarily receive suprathreshold input 
once per period that re-synchronizes possibly split- 
up oscillators. In general, a stable attractor has a 
contracting neighborhood in state space from which 
perturbed trajectories return to the attractor. If an 
attractor is unstable there are trajectories arbitrar- 
ily close to the attractor, which diverge from it, such 
that unstable attractors do not exhibit a contracting 
neighborhood. In other words, an unstable attractor 
has to be located at the boundary of its own basin. 
Therefore, the intuitive expectation is that parame- 
ters of the system have to be precisely specified in 
order to keep a periodic orbit simultaneously attract- 
ing and unstable. This leads to the question whether 
the physical parameters of the system, in particular 
£, r, and TV, need to be precisely tuned to obtain 
unstable attractors. 

To answer the question, how common unstable at- 
tractors actually are, we numerically estimated the 
fraction Pn{N) of state space occupied by basins of 
unstable attractors. To obtain this estimate, we ini- 
tialized the system with 1000 random initial phase 
vectors, drawn from the uniform distribution on 
[0,1]^. Whenever a period-one orbit was reached, 
we applied one random phase perturbation 8 drawn 
from the uniform distribution on [0,cr]^where we 
chose (T = 10~^, a value well below all scales that are 
determined by the model parameters, in particular 
a <^ e/N for network sizes up to iV « 10^. If per- 
turbed trajectories did not return to the original at- 
tractor, it was counted unstable. If no period-one or- 
bit was reached from a random initial state but, e.g., 
orbits of higher period, these were not tested for sta- 
bility. Thus, the numerical method used estimates 
a lower bound on Pu{N). As an example. Fig. [ll] a 
displays such an estimate of Pu{N) for e = 0.2 and 
T = 0.15. 

While for these parameters unstable attractors are 
absent if networks are too small (here iV < 4) and 
coexist with stable attractors in larger networks, the 
fraction approaches one for N ^ 1. Other param- 
eters (e — 0.2, T — 0.25) yield a different depen- 
dence on the network size N. Once again, unsta- 
ble attractors arise only if the network size is not 
too small {N > 3). Yet, the fraction PuiN) is only 
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Figure 11: Prevalence of unstable attractors in large net- 
works, (a) Unstable attractors prevail for large networks 
for certain parameters (e = 0.2, r = 0.15), but (b) are 
not important in large networks for other parameters 
(e = 0.2, r = 0.25). The fraction Pu{N) was estimated 
for every A'^ < 128 from 1000 random initial phase vec- 
tors, drawn from the uniform distribution on [0, l]'^. 




Delay r 

Figure 12: Unstable attractors persist in a wide region 
of parameter space in large networks (A'' = 100). Pa- 
rameters with Pu(lOO) > 0.5 are marked in black; here 
Pu(lOO) was estimated from 100 random initial phase vec- 
tors, drawn from the uniform distribution on [0, 1]^, for 
every set of parameters (resolution At = Ae = 0.005). 



substantial for moderate network sizes near iV « 10 
and approaches zero for large N. More generally, 
we observed that unstable attractors are absent in 
small networks (cf. ^ for the case = 2) and 
Pu{N) approaches either zero or one in large net- 
works, depending on the parameters. For networks 
of iV = 100 oscillators, Fig. ^2] shows the region of 
parameter space in which unstable attractors prevail 
(pu(lOO) > 0.5, estimated from 100 random initial 
phase vectors). As this region covers a substantial 
part of parameter space, precise parameter tuning is 
not needed to obtain unstable attractors. Further- 
more we find the same qualitative behavior indepen- 
dent of the detailed form of U(4>). This indicates that 
the occurrence of unstable attractors is a robust col- 
lective phenomenon in this model class of networks 
of excitatorily pulse-coupled oscillators. 

It is instructive to note that, on theoretical 
grounds, a period-one orbit is stable only if the clus- 
ters have a difference in firing times of At = r, such 
that its period is an integer multiple of r. The 
pulse sent by every single cluster then leads to a 
suprathreshold input to the following cluster, that 
in response to this input sends a pulse. If these con- 
ditions are not satisfied for all clusters of a period- 
one orbit, this orbit is unstable. In particular, it 
is unstable against a split-up of (at least) one clus- 
ter. The same orbit may also be attracting if such 
a cluster is formed in a region of state space that is 
located remote from the attractor as exemplified by 



the analysis in section VI 



Whereas the analysis presented above gives some 
insights into why unstable attractors exist and 
demonstrates that they prevail under variation of 
parameters, the precise reasons for their prevalence 
await discovery in future studies. 



VIII. CONCLUSIONS AND DISCUSSION 

The occurrence of unstable attractors per se is 
an intriguing phenomenon because it contradicts the 
common intuition about the stable nature of attract- 
ing invariant sets in dynamical systems. Our results 
suggest, that there are systems of pulse-coupled units 
in which unstable attractors may be the rule rather 
than the exception. 

Unstable attractors persist under various classes of 
structural modifications. For instance, preliminary 
studies on networks with randomly diluted connec- 
tivity suggest, that a symmetric, all-to-all connectiv- 
ity is not required ||3^. In addition, unstable attrac- 
tors also arise naturally in networks of inhibitorily 
coupled oscillators |Q, if a lower threshold is in- 
troduced and the function U is taken to be convex 
down, U"{4>) > 0, in a certain range of phase values. 
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a model variant motivated by experiments in certain 
biological neural systems | |29| |. 

Moreover, it is expected that every system ob- 
tained by a sufficiently small structural perturbation 
from the one considered here will exhibit a similar set 
of saddle periodic orbits, because linearly unstable 
states can generally not be stabilized by such a per- 
turbation. Although, in general, these orbits may no 
longer be attracting, their dynamical consequences 
are expected to persist. In particular, a switching 
along heterocHnic connections may occur in the pres- 
ence of noisy or deterministic, time-varying signals. 
As in the original system, the sequence of states 
reached may be determined by the directions into 
which such a signal guides the trajectory. By in- 
creasing and decreasing the strength of this signal, 
the time-scale of switching may be decreased and in- 
creased, respectively, due to the linear instability. 

Furthermore, switching among unstable states 
also occurs in systems of continuously phase-coupled 
oscillators js^, that can be obtained from pulse- 
coupled oscillators in a certain limit of weak cou- 
pHng JM. In particular. Hansel, Mato, and Meu- 
nier [|34|| show that a system of phase-coupled oscil- 
lators may switch back and forth among pairs of two- 
cluster states. Working in the limit of infinitely fast 
response, i.e. discontinuous phase jumps, we have 
demonstrated that far more compHcated switching 
transitions may occur in large networks if the oscil- 
lators are pulse-coupled. 

Unstable attractors add a high degree of fiexibility 
to a system allowing, e.g., switching from one attrac- 
tor towards a set of other attractors. This may be 
utilized for specific functions. If, for instance, the 
convergence of the state of the system towards an 



attractor has a functional role, such as the solution 
of a computational task ||l|, ^, ^, the flexibility in- 
duced by unstable attractors can provide the system 
with a unique advantage: If the attractor is stable, 
it will be hard to leave it after convergence, e.g. the 
completion of a task. With an unstable attractor, 
however, a small perturbation is sufficient to leave 
the attractor after convergence and proceed with the 
next task. This dynamical fiexibility might be used 
efficiently for the design of artificial systems and be 
highly advantageous to the computational capability 
of natural systems like neuronal networks. Interest- 
ingly, it has recently been shown that certain models 
of neural networks are capable of dynamically en- 
coding information as trajectories near heteroclinic 
connections ||3^ . 

In order to fully understand the capabilities of sys- 
tems exhibiting unstable attractors or unstable peri- 
odic orbits linked by heteroclinic connections and to 
learn to design such systems for specific functions it 
will be of major importance to further analyze the re- 
quirements for the occurrence of unstable attractors 
in dynamical systems and the factors which shape 
their basins. Our results indicate that networks of 
pulse-coupled oscillators may be a promising starting 
point for such investigations. 
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